﻿using System;
using System.Text;
using System.Drawing;
using System.Buffers;
using System.Collections;
using System.Collections.Generic;
using System.Runtime.InteropServices;

public static partial class NativeAOT
{
    [UnmanagedCallersOnly(EntryPoint = "simp")]
    public static unsafe double simp(double a, double b, double eps, IntPtr f_x_ptr)
    {
        f_x = Marshal.GetDelegateForFunctionPointer<delegatefunc_x>(f_x_ptr);

        return simp(a, b, eps);
    }

    /// <summary>
    /// Simpson求积法
    /// f计算被积函数值f(x)的函数名。
    /// </summary>
    /// <param name="a">积分下限。</param>
    /// <param name="b">积分上限。要求b>a。</param>
    /// <param name="eps">积分精度要求</param>
    /// <returns>函数返回积分值。</returns>
    public static unsafe double simp(double a, double b, double eps)
    {
        int n, k;
        double h, t1, t2, s1, s2 = 0.0, ep, p, x;

        n = 1;
        h = b - a;
        t1 = h * (f_x(a) + f_x(b)) / 2.0;
        s1 = t1;
        ep = eps + 1.0;
        while (ep >= eps)
        {
            p = 0.0;
            for (k = 0; k <= n - 1; k++)
            {
                x = a + (k + 0.5) * h;
                p = p + f_x(x);
            }
            t2 = (t1 + h * p) / 2.0;
            s2 = (4.0 * t2 - t1) / 3.0;
            ep = Math.Abs(s2 - s1);
            t1 = t2;
            s1 = s2;
            n = n + n;
            h = h / 2.0;
        }
        return (s2);
    }

    /*
      int main()
      { 
          double a,b,eps,t,simpf(double);
          a=0.0; b=1.0; eps=0.000001;
          t=simp(a,b,eps,simpf);
          cout <<"t = " <<t <<endl;
          return 0;
      }
    // 计算被积函数值
      double simpf(double x)
      { 
          double y;
          y=log(1.0+x)/(1.0+x*x);
          return(y);
      }
    */
}

